//------------------------------------------------------------------------------
// Fichiers d'ent�te utilises
//------------------------------------------------------------------------------ 
#include "OtsuThresholdFinder.hpp" 

namespace ImageAnalysis
{
	
	OtsuThresholdFinder::OtsuThresholdFinder():ThresholdFinder(){}
	
	OtsuThresholdFinder::~OtsuThresholdFinder(){}
	
	unsigned char * OtsuThresholdFinder::GetOneThreshold(Region* regionToStudy, Gradiant* gradiantToStudy)
	{
			
		int i,j,k=0;
			
		float sigma[256];
		float muT=0.0;
			
		Histogram * histo = gradiantToStudy->GetHistogram(regionToStudy);
	
		// calcul de muT
		muT = 0;
		for (i=0;i<256;i++)
			muT += (i+1)*histo->GetNumberOfPixel(i);
	
		muT /= (float)(regionToStudy->Width*regionToStudy->Height*1.0);
	
		// calcul des sigmas
		for (i=1; i<255; i++)
		{
			sigma[i] = ComputeSigmas(i,muT,regionToStudy->Width,regionToStudy->Height,histo);
		}
			
	
		// recherche du max
		i=0; 
		k=0;
		for (j=1; j<255;j++)
		{
			if (i+1 < sigma[j])   //on a un nouveau maximum
			{ 
				k=j; i=(int)sigma[j];  // on met le seuil dans k
			}
		}
	
		
		unsigned char * threshold = new unsigned char [1];
		threshold[0] = (unsigned char) k;
	
		delete histo;
		return threshold;
	}
	
	unsigned char * OtsuThresholdFinder::GetTwoThresholds(Region* regionToStudy, Gradiant* gradiantToStudy)
	{
		//on appelle la méthode pour un seuil pour trouver le seuil haut
		unsigned char * thresholdRight = GetOneThreshold( regionToStudy, gradiantToStudy);
	
		// Calcul de l'histogramme de la région
		Histogram * histo = gradiantToStudy->GetHistogram(regionToStudy);

		//on récupère le premier mode de l'histogramme
		//pour simplifier ici, on considère qu'il s'agit du premier pic de l'histogramme à contenir plus de 1/20ème des pixels
		int quantity = regionToStudy->Width*regionToStudy->Height/20;
		unsigned int modeIndex = 0;
		for(unsigned int i=0;i<=255;i++)
		{
			if(histo->GetNumberOfPixel(i)>=quantity)
			{
				modeIndex = i;
				break;
			}
		}
		
		//d'après l'article que nous avons trouvé, on obtient le seuil bas ainsi :
		unsigned char thresholdLeft = (unsigned char)(histo->GetNumberOfPixel(modeIndex) + 0.87f * (thresholdRight[0] - histo->GetNumberOfPixel(modeIndex)));
	
		unsigned char * thresholds = new unsigned char [2];
		thresholds[0] = thresholdLeft;
		thresholds[1] = *thresholdRight;
	
		delete thresholdRight;
		delete histo;
		return thresholds;
	}
	
	float OtsuThresholdFinder::ComputeSigmas(int k, float mu, int taille_x, int taille_y, Histogram* histo)
	{
		double tmp, w, tmp2;
		double N = (double)taille_x*taille_y;
		int i;
	
		w = 0;
		tmp = 0;
		for (i=0 ; i<=k; i++)
		{
			w += histo->GetNumberOfPixel(i);
			tmp += (i+1)*histo->GetNumberOfPixel(i);   //on calcule mu(k)
		}
		w /= N;
		tmp /= N;
		tmp2 = (mu*w-tmp)*(mu*w-tmp)*(w*(1-w));
	
		return (float)tmp2;
	}

}



